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PREFACE 


The  physical  assumption  of  laminar  flow  in  the  large 
blood  vessels  further  complicates  the  mathematical  model 
of  drug  distrioutlon  in  the  body  by  introducing  convolution 
terms  which  are  difficult  computationally.  This  Memorandum 
tackles  the  new  equations  of  the  chemotherapy  model  and 
presents  a  method  suitable  for  programming  this  model  as 
well  as  other  biological  systems  involving  equations  of 
this  type. 
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SUMMARY 

The  new  model  of  drug  distribution  In  the  body  In¬ 
corporates  the  exchange  between  the  stationary  and  flowing 
phases  in  the  large  blood  vessels.  This  introduces  com¬ 
putationally  difficult  convolution  terms.  The  method  of 
differential  approximation  applied  to  the  convolution 
equations  reduces  this  model  to  a  system  of  differential- 
difference  equations  which  can  be  solved  computationally. 
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A  NUMERICAL  APPROACH  TO  THE  CONVOLUTION  EQUATIONS  OF 
A  MATHEMATICAL  MODEL  OF  CHEMOTHERAPY 

1 .  INTRODUCTION 

Mathematical  models  of  drug  distribution  presented  in 
previous  papers  [l  —  5]  involved  the  physical  assumption  of 
volume  displacement  flow  in  the  large  blood  vessels  and  led 
to  a  system  of  differential— difference  equations  with  time 
delays  arising  from  the  non— zero  recirculation  time  of  the 
blood.  These  time  lags  together  with  the  parameters  in¬ 
volving  the  heart  as  a  mixing  pot  compensated  for  the 
assumption  of  "slug"  flow  with  no  mixing  in  the  arteries 
and  veins.  This  model  initiated  continuing  research  on 
the  numerical  solution  of  differential— difference  equations 
(see  [7,8,9])- 

The  work  of  Cooper  and  Jacquez  incorporates  the  exchange 
between  the  stationary  and  flowing  phases  in  the  large  blood 
vessels  using  some  results  of  Landahl  (see  [6],[l0]).  These 
assumptions  lead  to  a  set  of  differential  equations  together 
with  some  equations  containing  convolution  terms,  an  addi¬ 
tional  complication.  The  problem  of  obtaining  a  numerical 
solution  of  equations  of  this  nature  has  in  turn  stimulated 
considerable  research  on  various  approximation  techniques 
since  convolution  terms  introduce  significant  time  and 
storage  problems.  We  shall  apply  here  come  recent  results 
on  differential  approximation  (see  [11])  to  the  equations 
of  our  model  and  obtain  thereby  a  new  set  of  differential- 
difference  equations  which  can  be  solved  by  means  of  the 
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techniques  described  in  the  foregoing  references.  Further 
details  will  be  presented  in  a  subsequent  paper  describing 
the  digital  computer  program. 

2.  THE  EQUATIONS  OF  THE  MODEL 

The  equations  (see  [4],[l0])  of  our  two  organ  model 
consist  of  a  set  of  16  first  order  differential  equations 
for  the  functions  u.  .  v,,  w,,  z,,  j=l,  2,  m=2,  3,..., 6, 

J J  J  J 

and  the  following  3  convolution  equations  for  the  functions 
U1  l»  u2  1'  (un^ess  otherwise  specified  these  functions 
are  functions  of  t): 

t— s 

(1)  -  H1  +  d1  J  UL  ( q )  G  ( t— s— q )  dq 

o 

t — s 

(2)  ^2,1  =  H2  +  d£  J  UL(q)G(t-s-q)dq 

o 

t-s 

(3)  UL  =  H3  +  d3  J  u1^g(q)Cr(  t— s— q )  dq 

o 

t— s 

+  J  u2^(q)G(t— s— q)dq  . 

o 

here  s  =  ivf/c,  wl th  time  dimension, 

a  =  e  (-'-^/c)^  b  =  K2i/v  c,  h  =  k/v 

S  -  ^  s 

O-j^  —  ^ 1  ^  =  ^ 2  ^ 9  a3  =  ^lacl 9  =  ^lac2 9 
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nc ..  f  1  k  A  ,  f 

Sla  «L(t-a)  -  u1(1  -  -^-S-  Uljl  -  vx 


H£  -  — |  f2(t)  +  gja  UL(t-s) 


U2,l 


k  A  „ 
e  e2 


c  T_  „  ClUl,6(t-s)  +  c2U2.6^SK  _  ,  t t 

H3  =  V*  sla  - - - - -  +  g2fl(t)  ”  U 


a(x)  =  e  hX  )  =  e  hx  ,/b/x  I  (2,/bx) 

dx 


where  1^  is  the  1— st  modified  Bessel  function  of  the  first 
kind.  It  should  be  pointed  out  that  the  19  functions  in 
this  system  of  equations  are  zero  for  t  <  0. 

For  simplicity  the  model  assumes  the  same  length  l 
for  the  large  blood  vessels  leading  to  and  from  organs  1 
and  2.  This  leads  to  differential— difference  equations 
Involving  one  time  lag,  s.  We  can,  however,  handle  dif¬ 
ferent  commensurate  time  lags  (see  [12]).  This,  of  course, 
increases  time  and  storage  requirements. 

We  will  proceed  to  replace  each  of  the  three  convolu¬ 
tion  equations  by  a  set  of  differential— difference  equations 
using  a  procedure  described  in  the  sections  that  follow. 

3.  DIFFERENTIAL  APPROXIMATION 

The  method  of  differential  approximation  (see  [ll]) 
permits  us  to  obtain  a  set  of  coefficients,  a.,  i  =  1,...,N, 

\  T 

so  that  G^“^(x)  +  Z  a.G^J— ^(x)  =  0,  in  the  mean  square 

i=l  1 

sense.  We  proceed  a.s  follows. 


The  modified  Bessel  function  of  the  first  kind  satisfies 


the  second  order  differential  equation 

I^(z)  +  \  I^(z)  —  (1  +  -i>)  I^(z)  =  0  . 

Prom  this  it  follows  that 

xg"(x)  +  G*  (x) (2hx+2 )  +  G(x) (h2x+2h-b)  =  0  . 

By  repeated  differentiation  we  obtain 

x(j(N+1)  +  (2hx+N+l)G^N^  +  ( 2hN+h2x— b )  G  ^  1  ^ 

+  h2(N-l)G(N~2)  =  0  . 

Let  G^Cx)  =  G^^(x)  and  consider  the  N+l  differential 
equations 


3x  Gi+1 


1  =  0,  1,  .  .  N— 1 


^  |  |(2hx+N+l)GN  +  (2hIJ-(h2x-b)G1^1  +  h2  (N-ljG^g} 

The  initial  conditions  G^(0)  can  be  obtained  from  the 
known  series  expansion  for  I-^(x),  namely 

-W_  .  b2x  .  b3x2  .  .  b<n+1>xn  .  ^ 

G(x)=e  1  b  +  —rr—  +  - —  +.  .  .  +  - 5  +.  .  .  J . 

V  2  3(21)2  (n+l)(n!)2  7 


so  that  over 


We  wish  to  determine  the  coefficients 
a  suitable  range  T  v/e  obtain  the  minimum 


=  min  /  [G„(x)  +  iSia1G(N_l)J  dx  . 


Setting  the  partial  derivatives  of  M  with  respect  to  a^ 
equal  to  zero,  we  obtain  the  linear  system 


N  1  1 

Z  a.  G..  .G„  ,  dx  =  - 

<=1  J  J  N-i  J 


GM-iaN  ’ 


i  =  1,  N 


-L 


k,  q  =  0,  . . . ,  N  —  1 


J. 

Sc  =  /  GkGNdx  > 


k  =  0,  .  .  . ,  N  -  1 


To  the  set  of  equations  (4),  add  the  following  differential 
equations 


k,q  _ 


=  ^  , 


=  ~  GkGM  ' 


q ( 0)  —  0,  k, q  —  0,,..., 


0^(0)  =0,  k  =  0,  .  N  - 


Integrating  the  combined  set  of  equations  (4)  and  (5) 
and  taking  the  final  values  of  q  and  at  x  ■  T,  we 
solve  the  following  system  of  N  linear  equations  in  N 
unknowns  aN_q»  q  =  0,  . . N  —  1, 

N— 3. 

q=0  aN-qPk,q^T^  "  •  k  "  °>  N  ~ 

Ihe  solution  gives  us  a^,  ai  • 

4.  REDUCTION  TO  A  SYSTEM  OF  DIFFERENTIAL-DIFFERENCE 

— Tggpgrc ms - 

For  simplicity  of  presentation,  let  us  assume  N  =  3 
for  the  order  of  the  differential  approximation.  Let  us 
define 

t 

X(t)  =■  j  UL(q)G(t-q)dq 
o 

t 

Y(t)  =  J  U1#g(q)G(t-q)dq 
o 

t 

Z(t)  =  J  U^6(q)G(t-q)dq  , 
o 

and 


Ihe  3  convolution  equations  become 
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(6) 

u-j-  =  "*■  d^X(t— s) 

(7) 

u2  i  =  H2  +  dgX(t-s) 

(S) 

1 

UL  =  H3  +  d3Y(t~s)  +  d4Z(t-s) 

To 

these  equations  adjoin  the  following 

set : 

r 

x'  =  X]_, 

X(0)  =  0 

x|  -  x2. 

xx(o)  =  0 

X2  =  h4(x,x1,x2,ul,u[,u^). 

x2(o)  =  0 

y'  -  V 

Y  (  0 )  =  0 

(9) 

< 

Yi  -  V 

Yx(0)  -  0 

'  /  1  " 

Y2  =  ( Y,  Y-^ ,  Y2  ,  u.^  t  g ,  u ^  g ,  u  1  g - ) , 

y2(o)  =  0 

z'  =  z1. 

Z(C)  =  0 

1 

Z1  =  Z2 ’ 

Z1(c)  =  0 

Z2  =  iX4  (Z,2I,22,U2,C,U2,6,U2,6^  * 

Z2(0)  =  0 

Title 

function  and  the  Initial  conditions  ar 

e  obtained 

as 

follows:  Illustrating  with  X(t),  by  using 

Le i bn iz's 

rule  repeatedly,  we  get 

t 

(10)  x'(t)  =  a(0)UL(t)  +  J  g' (t-q)UL(q)dq 

o 

t 

(11)  X  (z)  =  G(0)U^(t)  +  a'(0)dL(t)  +  J  a"(t-q)UL(q)dq 

o 
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(12)  x"'(t)  =  G(0)u^(t)  +  g' (o)u^(t)  +  a"(o)uL(t) 

t 

+  y  G  '  (t-q)UL(q)dq  . 
o 

Combining  these  together  with  the  3  coefficients  a±  of 
our  differential  approximation,  we  obtain 


X  '  +a-LX  +a2x'+a3X  =  G(0)U^(t)+G' (OjU^tJ+G '(0)UL(t) 

+  a^ ( G ( 0 ) ( t ) +G ' ( 0 ) UL ( t ) )  +  a2G(0)UL(t) 

+  U^(q)(G  (t— q)  +  a^G  (t— q)+a2G  (t— q)+a3G(t— q)  )dq 


By  a  change  of  variable  and  the  result  of  Sec.  3,  the 
Integrand  in  the  last  term  vanishes.  The  function 
by  virtue  of  (13)  is  given  by 


( X,X-^,  X2  ,  Uj^,  —  — a3^— a2^1 — a1^2 


-|g  (0)+a1Gl  (0)+a2G(0)JuIj  +  g' (0)+a1G(0)  ju^  +  G(0)U: 


By  going  back  to  the  equations  of  the  model  (see  [4], [10]) 

i  ii  i  ii 

the  derivative  terms  UL>ui/uj  6,uj  in  the  function 

H^,  can  be  expressed  in  terms  of  the  original  functions  of 
our  set  of  28  differential-difference  equations,  for  a 


_o_ 


tv/o  organ  model  with  a  third  order  differential  approxi¬ 
mation  for  Gr(x).  Differential  approximation  of  order  k 
involves  a  system  of  19  +  3k  equations. 
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